pacman::p_load(sf,tmap,spdep,sf,tidyverse, dplyr, mapview)Take_home_Ex1
Overview
Objectives
Load Packages and Data
Load packages
Load data
Import aspiatial data
odbus <-read_csv("data/aspatial/origin_destination_bus_202310.csv")Getting the data of the different duration
weekdayAM_6_9<-odbus %>%
filter(DAY_TYPE == "WEEKDAY") %>%
filter(TIME_PER_HOUR >= 6 &
TIME_PER_HOUR <= 9) %>%
group_by(ORIGIN_PT_CODE) %>%
summarise(TRIPS = sum(TOTAL_TRIPS))
weekdayPM_5_8<-odbus %>%
filter(DAY_TYPE == "WEEKDAY") %>%
filter(TIME_PER_HOUR >= 17 &
TIME_PER_HOUR <= 20) %>%
group_by(ORIGIN_PT_CODE) %>%
summarise(TRIPS = sum(TOTAL_TRIPS))
weekendAM_11_14<-odbus %>%
filter(DAY_TYPE == "WEEKENDS/HOLIDAY") %>%
filter(TIME_PER_HOUR >= 11 &
TIME_PER_HOUR <= 14) %>%
group_by(ORIGIN_PT_CODE) %>%
summarise(TRIPS = sum(TOTAL_TRIPS))
weekendPM_4_7<-odbus %>%
filter(DAY_TYPE == "WEEKENDS/HOLIDAY") %>%
filter(TIME_PER_HOUR >= 16 &
TIME_PER_HOUR <= 19) %>%
group_by(ORIGIN_PT_CODE) %>%
summarise(TRIPS = sum(TOTAL_TRIPS))Import geospatial
busstop <- st_read(dsn = "data/geospatial",
layer = "BusStop") %>%
st_transform(crs = 3414)Reading layer `BusStop' from data source
`C:\Users\Lian Khye\Desktop\MITB\Geospatial\geartooth\ISSS624\Take_home_Ex1\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 5161 features and 3 fields
Geometry type: POINT
Dimension: XY
Bounding box: xmin: 3970.122 ymin: 26482.1 xmax: 48284.56 ymax: 52983.82
Projected CRS: SVY21glimpse(busstop)Rows: 5,161
Columns: 4
$ BUS_STOP_N <chr> "22069", "32071", "44331", "96081", "11561", "66191", "2338…
$ BUS_ROOF_N <chr> "B06", "B23", "B01", "B05", "B05", "B03", "B02A", "B02", "B…
$ LOC_DESC <chr> "OPP CEVA LOGISTICS", "AFT TRACK 13", "BLK 239", "GRACE IND…
$ geometry <POINT [m]> POINT (13576.31 32883.65), POINT (13228.59 44206.38),…mpsz <- st_read(dsn = "data/geospatial",
layer = "MPSZ-2019") %>%
st_transform(crs = 3414)Reading layer `MPSZ-2019' from data source
`C:\Users\Lian Khye\Desktop\MITB\Geospatial\geartooth\ISSS624\Take_home_Ex1\data\geospatial'
using driver `ESRI Shapefile'
Simple feature collection with 332 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 103.6057 ymin: 1.158699 xmax: 104.0885 ymax: 1.470775
Geodetic CRS: WGS 84busstop_mpsz <- st_intersection(busstop, mpsz) %>%
select(BUS_STOP_N, SUBZONE_C) %>%
st_drop_geometry()AM_6_9 <- left_join(weekdayAM_6_9 , busstop_mpsz,
by = c("ORIGIN_PT_CODE" = "BUS_STOP_N")) %>%
rename(ORIGIN_BS = ORIGIN_PT_CODE,
ORIGIN_SZ = SUBZONE_C)
PM_5_8 <- left_join(weekdayPM_5_8 , busstop_mpsz,
by = c("ORIGIN_PT_CODE" = "BUS_STOP_N")) %>%
rename(ORIGIN_BS = ORIGIN_PT_CODE,
ORIGIN_SZ = SUBZONE_C)
AM_11_2<- left_join(weekendAM_11_14 , busstop_mpsz,
by = c("ORIGIN_PT_CODE" = "BUS_STOP_N")) %>%
rename(ORIGIN_BS = ORIGIN_PT_CODE,
ORIGIN_SZ = SUBZONE_C)
PM_4_7<- left_join(weekendPM_4_7 , busstop_mpsz,
by = c("ORIGIN_PT_CODE" = "BUS_STOP_N")) %>%
rename(ORIGIN_BS = ORIGIN_PT_CODE,
ORIGIN_SZ = SUBZONE_C)duplicate1 <- AM_6_9 %>%
group_by_all() %>%
filter(n()>1) %>%
ungroup()
duplicate2 <- PM_5_8 %>%
group_by_all() %>%
filter(n()>1) %>%
ungroup()
duplicate3 <- AM_11_2 %>%
group_by_all() %>%
filter(n()>1) %>%
ungroup()
duplicate4 <- PM_5_8 %>%
group_by_all() %>%
filter(n()>1) %>%
ungroup()c(duplicate1,duplicate2,duplicate3,duplicate4)$ORIGIN_BS
[1] "11009" "11009" "22501" "22501" "43709" "43709" "47201" "47201" "51071"
[10] "51071" "53041" "53041" "62251" "62251" "67421" "67421" "68091" "68091"
[19] "77329" "77329" "82221" "82221" "96319" "96319" "97079" "97079"
$TRIPS
[1] 17311 17311 11384 11384 1329 1329 25910 25910 6698 6698 4462 4462
[13] 8605 8605 258 258 254 254 337 337 3998 3998 675 675
[25] 3822 3822
$ORIGIN_SZ
[1] "QTSZ01" "QTSZ01" "JWSZ09" "JWSZ09" "BKSZ07" "BKSZ07" "WDSZ07" "WDSZ07"
[9] "CCSZ01" "CCSZ01" "BSSZ01" "BSSZ01" "HGSZ03" "HGSZ03" "SESZ04" "SESZ04"
[17] "SLSZ04" "SLSZ04" "PRSZ03" "PRSZ03" "GLSZ05" "GLSZ05" "TMSZ05" "TMSZ05"
[25] "CHSZ02" "CHSZ02"
$ORIGIN_BS
[1] "11009" "11009" "22501" "22501" "43709" "43709" "47201" "47201" "51071"
[10] "51071" "53041" "53041" "62251" "62251" "67421" "67421" "68091" "68091"
[19] "77329" "77329" "82221" "82221" "96319" "96319" "97079" "97079"
$TRIPS
[1] 6229 6229 2533 2533 1608 1608 15689 15689 14145 14145 10621 10621
[13] 1874 1874 297 297 982 982 326 326 2871 2871 2996 2996
[25] 1980 1980
$ORIGIN_SZ
[1] "QTSZ01" "QTSZ01" "JWSZ09" "JWSZ09" "BKSZ07" "BKSZ07" "WDSZ07" "WDSZ07"
[9] "CCSZ01" "CCSZ01" "BSSZ01" "BSSZ01" "HGSZ03" "HGSZ03" "SESZ04" "SESZ04"
[17] "SLSZ04" "SLSZ04" "PRSZ03" "PRSZ03" "GLSZ05" "GLSZ05" "TMSZ05" "TMSZ05"
[25] "CHSZ02" "CHSZ02"
$ORIGIN_BS
[1] "11009" "11009" "22501" "22501" "43709" "43709" "47201" "47201" "51071"
[10] "51071" "53041" "53041" "62251" "62251" "67421" "67421" "68091" "68091"
[19] "77329" "77329" "82221" "82221" "96319" "96319" "97079" "97079"
$TRIPS
[1] 4698 4698 1539 1539 415 415 2204 2204 5163 5163 4001 4001 1431 1431 91
[16] 91 223 223 92 92 2249 2249 164 164 290 290
$ORIGIN_SZ
[1] "QTSZ01" "QTSZ01" "JWSZ09" "JWSZ09" "BKSZ07" "BKSZ07" "WDSZ07" "WDSZ07"
[9] "CCSZ01" "CCSZ01" "BSSZ01" "BSSZ01" "HGSZ03" "HGSZ03" "SESZ04" "SESZ04"
[17] "SLSZ04" "SLSZ04" "PRSZ03" "PRSZ03" "GLSZ05" "GLSZ05" "TMSZ05" "TMSZ05"
[25] "CHSZ02" "CHSZ02"
$ORIGIN_BS
[1] "11009" "11009" "22501" "22501" "43709" "43709" "47201" "47201" "51071"
[10] "51071" "53041" "53041" "62251" "62251" "67421" "67421" "68091" "68091"
[19] "77329" "77329" "82221" "82221" "96319" "96319" "97079" "97079"
$TRIPS
[1] 6229 6229 2533 2533 1608 1608 15689 15689 14145 14145 10621 10621
[13] 1874 1874 297 297 982 982 326 326 2871 2871 2996 2996
[25] 1980 1980
$ORIGIN_SZ
[1] "QTSZ01" "QTSZ01" "JWSZ09" "JWSZ09" "BKSZ07" "BKSZ07" "WDSZ07" "WDSZ07"
[9] "CCSZ01" "CCSZ01" "BSSZ01" "BSSZ01" "HGSZ03" "HGSZ03" "SESZ04" "SESZ04"
[17] "SLSZ04" "SLSZ04" "PRSZ03" "PRSZ03" "GLSZ05" "GLSZ05" "TMSZ05" "TMSZ05"
[25] "CHSZ02" "CHSZ02"AM_6_9 <- unique(AM_6_9)
PM_5_8 <- unique(PM_5_8)
AM_11_2 <- unique(AM_11_2)
PM_4_7 <- unique(PM_4_7)Check if duplicates cleared
mpsz_AM_6_9 <- left_join(mpsz,AM_6_9,by = c("SUBZONE_C" = "ORIGIN_SZ"))
mpsz_PM_5_8 <- left_join(mpsz,PM_5_8,by = c("SUBZONE_C" = "ORIGIN_SZ"))
mpsz_AM_11_2 <- left_join(mpsz,AM_11_2,by = c("SUBZONE_C" = "ORIGIN_SZ"))
mpsz_PM_4_7 <- left_join(mpsz,PM_4_7 ,by = c("SUBZONE_C" = "ORIGIN_SZ"))tmap_mode("plot")
AM_6_9_map <- tm_shape(mpsz_AM_6_9)+
tm_fill("TRIPS",
style = "quantile",
palette = "Blues",
title = "Passenger trips") +
tm_layout(main.title = "Passenger trips generated Weekday 6am-9am",
main.title.position = "center",
main.title.size = 0.7,
legend.height = 0.45,
legend.width = 0.35,
frame = TRUE) +
tm_borders(alpha = 0.5) +
tm_compass(type="8star", size = 1) +
tm_scale_bar() +
tm_grid(alpha =0.2) +
tm_credits("Source: Planning Sub-zone boundary from URA\n and Passenger trips data from LTA",
position = c("left", "bottom"))
PM_5_8_map <- tm_shape(mpsz_PM_5_8)+
tm_fill("TRIPS",
style = "quantile",
palette = "Greens",
title = "Passenger trips") +
tm_layout(main.title = "Passenger trips generated Weekday 5pm-8pm",
main.title.position = "center",
main.title.size = 0.7,
legend.height = 0.45,
legend.width = 0.35,
frame = TRUE) +
tm_borders(alpha = 0.5) +
tm_compass(type="8star", size = 1) +
tm_scale_bar() +
tm_grid(alpha =0.2) +
tm_credits("Source: Planning Sub-zone boundary from URA\n and Passenger trips data from LTA",
position = c("left", "bottom"))
AM_11_2_map <- tm_shape(mpsz_AM_11_2)+
tm_fill("TRIPS",
style = "quantile",
palette = "Reds",
title = "Passenger trips") +
tm_layout(main.title = "Passenger trips generated Weekend/holidays 11am-2pm",
main.title.position = "center",
main.title.size = 0.7,
legend.height = 0.45,
legend.width = 0.35,
frame = TRUE) +
tm_borders(alpha = 0.5) +
tm_compass(type="8star", size = 1) +
tm_scale_bar() +
tm_grid(alpha =0.2) +
tm_credits("Source: Planning Sub-zone boundary from URA\n and Passenger trips data from LTA",
position = c("left", "bottom"))
PM_4_7_map <- tm_shape(mpsz_PM_4_7)+
tm_fill("TRIPS",
style = "quantile",
palette = "Greens",
title = "Passenger trips") +
tm_layout(main.title = "Passenger trips generated Weekends/holidays 4pm-7pm",
main.title.position = "center",
main.title.size = 0.7,
legend.height = 0.45,
legend.width = 0.35,
frame = TRUE) +
tm_borders(alpha = 0.5) +
tm_compass(type="8star", size = 1) +
tm_scale_bar() +
tm_grid(alpha =0.2) +
tm_credits("Source: Planning Sub-zone boundary from URA\n and Passenger trips data from LTA",
position = c("left", "bottom"))#tmap_arrange(AM_6_9_map, PM_5_8_map,AM_11_2_map,PM_4_7_map, asp = 2, ncol = 2)Transform data to hex friendly
gtest_points = mpsz_AM_6_9 %>%
# lng/lat value are missing in some records
filter(!is.na(geometry)) %>%
st_as_sf(coords = "geometry", crs = 3414, remove = FALSE)Making grids
hex_grid <- st_make_grid(gtest_points, c(750, 750), what = "polygons", square = FALSE)
grid_sf = st_sf(hex_grid) %>%
# add grid ID
mutate(grid_id = 1:length(lengths(hex_grid)))contents <- st_intersects(hex_grid, gtest_points)
for (i in 1:length(contents)) {
x = 0
for (j in unlist(i)){
x = x + gtest_points$TRIPS[j]
}
grid_sf$ridership[i] = x
}grid_sf$ridership = lengths(st_intersects(hex_grid, gtest_points))
grid_count = filter(grid_sf, ridership > 0)AM_6_9_map <- tm_shape(grid_count)+
tm_fill("ridership",
style = "quantile",
palette = "Blues",
title = "Passenger trips") +
tm_layout(main.title = "Passenger trips generated Weekday 6am-9am",
main.title.position = "center",
main.title.size = 0.7,
legend.height = 0.45,
legend.width = 0.35,
frame = TRUE) +
tm_borders(alpha = 0.5) +
tm_compass(type="8star", size = 1) +
tm_scale_bar() +
tm_grid(alpha =0.2) +
tm_credits("Source: Planning Sub-zone boundary from URA\n and Passenger trips data from LTA",
position = c("left", "bottom"))
AM_6_9_map
tmap_mode("view")
map_honeycomb = tm_shape(grid_count) +
tm_fill(
col = "ridership",
palette = "Blues",
style = "cont",
title = "Number of bus stops",
id = "grid_id",
showNA = FALSE,
alpha = 0.6,
popup.vars = c(
"Number of bus stops: " = "ridership"
),
popup.format = list(
ridership = list(format = "f", digits = 0)
)
) +
tm_borders(col = "grey40", lwd = 0.7)
map_honeycombmode(map_honeycomb)[1] "list"mode(map_honeycomb)[1] "list"mapview_test_points = mapview(gtest_points, cex = 3, alpha = .5, popup = NULL)
mapview_test_points